##########################################################
# Bootstrap functions
##########################################################
library(boot)

pabound.boot.lsdv <- function(data, index, x, y){
    X <- data[index, x]
    Y <- data[index, y]
    weapon <- as.factor(data[index, "weapon"])
    numseratt <- as.factor(data[index, "numseratt"])
    d <- na.omit(data.frame(X,Y,weapon,numseratt))
    if (length(unique(X))==1){
        res.pearl <- NA; res.exog <- c(NA,NA)
    } else {
        if(nlevels(d$weapon)==1 & nlevels(d$numseratt)>1) f <- Y ~ X + numseratt
        else if (nlevels(d$weapon)>1 & nlevels(d$numseratt)==1) f <- Y ~ X + weapon
        else if (nlevels(d$weapon)==1 & nlevels(d$numseratt)==1) f <- Y ~ X
        else f <- Y ~ X + weapon + numseratt
        DD <- data.frame(Y = d$Y, model.matrix(f)[,-1])
        fDD <- as.formula(paste(names(DD)[1], "~", paste(names(DD)[-1], collapse="+")))
        fit <- lm(fDD, DD)
        X0 <- X1 <- model.matrix(fit)
        X1[,2] <- 1
        X0[,2] <- 0
        P1 <- X1 %*% coef(fit)
        P0 <- X0 %*% coef(fit)
        P1 <- ifelse(P1 < 0 | P1 > 1, NA, P1)
        P0 <- ifelse(P0 < 0 | P0 > 1, NA, P0)
        res.pearl <- mean(pmax(0, P1-P0)/P1, na.rm=T)
        res.exog <- c(mean(pmax(0, P1-P0)/P1, na.rm=T), 
            mean(pmin(P1, 1-P0)/P1, na.rm=T))
    }
    c(res.pearl, res.exog)
}

##########################################################
# Read data
##########################################################
load("olken-data.RData")

# Corrections to (presumed) original coding error
D$numseratt <- ifelse(D$numseratt >= 3, 3, D$numseratt)
D[D$cowcode=="YAR" & D$year==1978,]$numseratt <- 1
D[D$cowcode=="YUG" & D$year==1934,]$numseratt <- 1

# Remove non-serious attempts
D <- subset(D, seriousattempt == 1)

###########################################################
# Analysis of Regime Change
###########################################################

# Subsetting data
D.dem <- subset(D, democracy.l == 1)
D.aut <- subset(D, democracy.l == 0)
D.old <- subset(D, newleader == 0)
D.new <- subset(D, newleader == 1)
D.aut.new <- subset(D, democracy.l == 0 & newleader == 1)
D.aut.old <- subset(D, democracy.l == 0 & newleader == 0)

# Assassination and Regime Change
resdem.all.lsdv <- boot(D, pabound.boot.lsdv, 2000, x="success", y="democracy.absch")
resdem.AtoD.lsdv <- boot(D.aut, pabound.boot.lsdv, 2000, x="success", y="democracy.absch")
resdem.DtoA.lsdv <- boot(D.dem, pabound.boot.lsdv, 2000, x="success", y="democracy.absch")

# Duration by leader tenure
# All
resdem.AtoD.allf10.lsdv <- boot(D.aut, pabound.boot.lsdv, 2000, x="success", y="democracy.absch10")
resdem.AtoD.allf20.lsdv <- boot(D.aut, pabound.boot.lsdv, 2000, x="success", y="democracy.absch20")

# <= 10
resdem.AtoD.new.lsdv <- boot(D.aut.new, pabound.boot.lsdv, 2000, x="success", y="democracy.absch")
resdem.AtoD.newf10.lsdv <- boot(D.aut.new, pabound.boot.lsdv, 2000, x="success", y="democracy.absch10")
resdem.AtoD.newf20.lsdv <- boot(D.aut.new, pabound.boot.lsdv, 2000, x="success", y="democracy.absch20")

# > 10
resdem.AtoD.old.lsdv <- boot(D.aut.old, pabound.boot.lsdv, 2000, x="success", y="democracy.absch")
resdem.AtoD.oldf10.lsdv <- boot(D.aut.old, pabound.boot.lsdv, 2000, x="success", y="democracy.absch10")
resdem.AtoD.oldf20.lsdv <- boot(D.aut.old, pabound.boot.lsdv, 2000, x="success", y="democracy.absch20")

###################################################
# Discard inconsistent draws if any
###################################################
resdem.all.lsdv <- na.omit(resdem.all.lsdv$t)
resdem.AtoD.lsdv <- na.omit(resdem.AtoD.lsdv$t)
resdem.DtoA.lsdv <- na.omit(resdem.DtoA.lsdv$t)
resdem.AtoD.allf10.lsdv <- na.omit(resdem.AtoD.allf10.lsdv$t)
resdem.AtoD.allf20.lsdv <- na.omit(resdem.AtoD.allf20.lsdv$t)
resdem.AtoD.new.lsdv <- na.omit(resdem.AtoD.new.lsdv$t)
resdem.AtoD.newf10.lsdv <- na.omit(resdem.AtoD.newf10.lsdv$t)
resdem.AtoD.newf20.lsdv <- na.omit(resdem.AtoD.newf20.lsdv$t)
resdem.AtoD.old.lsdv <- na.omit(resdem.AtoD.old.lsdv$t)
resdem.AtoD.oldf10.lsdv <- na.omit(resdem.AtoD.oldf10.lsdv$t)
resdem.AtoD.oldf20.lsdv <- na.omit(resdem.AtoD.oldf20.lsdv$t)

###################################################
# Point estimates
###################################################
resdem.all.lsdv.pearl <- median(resdem.all.lsdv[,1])
resdem.AtoD.lsdv.pearl <- median(resdem.AtoD.lsdv[,1])
resdem.DtoA.lsdv.pearl <- median(resdem.DtoA.lsdv[,1])
resdem.AtoD.allf10.lsdv.pearl <- median(resdem.AtoD.allf10.lsdv[,1])
resdem.AtoD.allf20.lsdv.pearl <- median(resdem.AtoD.allf20.lsdv[,1])
resdem.AtoD.new.lsdv.pearl <- median(resdem.AtoD.new.lsdv[,1])
resdem.AtoD.newf10.lsdv.pearl <- median(resdem.AtoD.newf10.lsdv[,1])
resdem.AtoD.newf20.lsdv.pearl <- median(resdem.AtoD.newf20.lsdv[,1])
resdem.AtoD.old.lsdv.pearl <- median(resdem.AtoD.old.lsdv[,1])
resdem.AtoD.oldf10.lsdv.pearl <- median(resdem.AtoD.oldf10.lsdv[,1])
resdem.AtoD.oldf20.lsdv.pearl <- median(resdem.AtoD.oldf20.lsdv[,1])

resdem.all.lsdv.exog <- apply(resdem.all.lsdv[,2:3], 2, median)
resdem.AtoD.lsdv.exog <- apply(resdem.AtoD.lsdv[,2:3], 2, median)
resdem.DtoA.lsdv.exog <- apply(resdem.DtoA.lsdv[,2:3], 2, median)
resdem.AtoD.allf10.lsdv.exog <- apply(resdem.AtoD.allf10.lsdv[,2:3], 2, median)
resdem.AtoD.allf20.lsdv.exog <- apply(resdem.AtoD.allf20.lsdv[,2:3], 2, median)
resdem.AtoD.new.lsdv.exog <- apply(resdem.AtoD.new.lsdv[,2:3], 2, median)
resdem.AtoD.newf10.lsdv.exog <- apply(resdem.AtoD.newf10.lsdv[,2:3], 2, median)
resdem.AtoD.newf20.lsdv.exog <- apply(resdem.AtoD.newf20.lsdv[,2:3], 2, median)
resdem.AtoD.old.lsdv.exog <- apply(resdem.AtoD.old.lsdv[,2:3], 2, median)
resdem.AtoD.oldf10.lsdv.exog <- apply(resdem.AtoD.oldf10.lsdv[,2:3], 2, median)
resdem.AtoD.oldf20.lsdv.exog <- apply(resdem.AtoD.oldf20.lsdv[,2:3], 2, median)

################################################################################
# Confidence Intervals
################################################################################
source("IMCI.R")

CI.all.lsdv.pearl <- c(resdem.all.lsdv.pearl - qnorm(.975)*sd(resdem.all.lsdv[,1]),
    resdem.all.lsdv.pearl + qnorm(.975)*sd(resdem.all.lsdv[,1]))
CI.AtoD.lsdv.pearl <- c(resdem.AtoD.lsdv.pearl - qnorm(.975)*sd(resdem.AtoD.lsdv[,1]),
    resdem.AtoD.lsdv.pearl + qnorm(.975)*sd(resdem.AtoD.lsdv[,1]))
CI.DtoA.lsdv.pearl <- c(resdem.DtoA.lsdv.pearl - qnorm(.975)*sd(resdem.DtoA.lsdv[,1]),
    resdem.DtoA.lsdv.pearl + qnorm(.975)*sd(resdem.DtoA.lsdv[,1]))

CI.all.lsdv.pearl <- pmax(0, pmin(1, CI.all.lsdv.pearl))
CI.AtoD.lsdv.pearl <- pmax(0, pmin(1, CI.AtoD.lsdv.pearl))
CI.DtoA.lsdv.pearl <- pmax(0, pmin(1, CI.DtoA.lsdv.pearl))

CI.all.lsdv.exog <- IMCI(l=resdem.all.lsdv.exog[1], u=resdem.all.lsdv.exog[2],
    var.l=var(resdem.all.lsdv[,2], na.rm=T), var.u=var(resdem.all.lsdv[,3], na.rm=T))
CI.AtoD.lsdv.exog <- IMCI(l=resdem.AtoD.lsdv.exog[1], u=resdem.AtoD.lsdv.exog[2],
    var.l=var(resdem.AtoD.lsdv[,2], na.rm=T), var.u=var(resdem.AtoD.lsdv[,3], na.rm=T))
CI.DtoA.lsdv.exog <- IMCI(l=resdem.DtoA.lsdv.exog[1], u=resdem.DtoA.lsdv.exog[2],
    var.l=var(resdem.DtoA.lsdv[,2], na.rm=T), var.u=var(resdem.DtoA.lsdv[,3], na.rm=T))

CI.all.lsdv.exog <- pmax(0, pmin(1, CI.all.lsdv.exog$ci))
CI.AtoD.lsdv.exog <- pmax(0, pmin(1, CI.AtoD.lsdv.exog$ci))
CI.DtoA.lsdv.exog <- pmax(0, pmin(1, CI.DtoA.lsdv.exog$ci))

CI.AtoD.allf10.lsdv.pearl <- c(resdem.AtoD.allf10.lsdv.pearl - qnorm(.975)*sd(resdem.AtoD.allf10.lsdv[,1]),
    resdem.AtoD.allf10.lsdv.pearl + qnorm(.975)*sd(resdem.AtoD.allf10.lsdv[,1]))
CI.AtoD.allf20.lsdv.pearl <- c(resdem.AtoD.allf20.lsdv.pearl - qnorm(.975)*sd(resdem.AtoD.allf20.lsdv[,1]),
    resdem.AtoD.allf20.lsdv.pearl + qnorm(.975)*sd(resdem.AtoD.allf20.lsdv[,1]))
CI.AtoD.new.lsdv.pearl <- c(resdem.AtoD.new.lsdv.pearl - qnorm(.975)*sd(resdem.AtoD.new.lsdv[,1]),
    resdem.AtoD.new.lsdv.pearl + qnorm(.975)*sd(resdem.AtoD.new.lsdv[,1]))
CI.AtoD.newf10.lsdv.pearl <- c(resdem.AtoD.newf10.lsdv.pearl - qnorm(.975)*sd(resdem.AtoD.newf10.lsdv[,1]),
    resdem.AtoD.newf10.lsdv.pearl + qnorm(.975)*sd(resdem.AtoD.newf10.lsdv[,1]))
CI.AtoD.newf20.lsdv.pearl <- c(resdem.AtoD.newf20.lsdv.pearl - qnorm(.975)*sd(resdem.AtoD.newf20.lsdv[,1]),
    resdem.AtoD.newf20.lsdv.pearl + qnorm(.975)*sd(resdem.AtoD.newf20.lsdv[,1]))
CI.AtoD.old.lsdv.pearl <- c(resdem.AtoD.old.lsdv.pearl - qnorm(.975)*sd(resdem.AtoD.old.lsdv[,1]),
    resdem.AtoD.old.lsdv.pearl + qnorm(.975)*sd(resdem.AtoD.old.lsdv[,1]))
CI.AtoD.oldf10.lsdv.pearl <- c(resdem.AtoD.oldf10.lsdv.pearl - qnorm(.975)*sd(resdem.AtoD.oldf10.lsdv[,1]),
    resdem.AtoD.oldf10.lsdv.pearl + qnorm(.975)*sd(resdem.AtoD.oldf10.lsdv[,1]))
CI.AtoD.oldf20.lsdv.pearl <- c(resdem.AtoD.oldf20.lsdv.pearl - qnorm(.975)*sd(resdem.AtoD.oldf20.lsdv[,1]),
    resdem.AtoD.oldf20.lsdv.pearl + qnorm(.975)*sd(resdem.AtoD.oldf20.lsdv[,1]))

CI.AtoD.allf10.lsdv.pearl <- pmax(0, pmin(1, CI.AtoD.allf10.lsdv.pearl))
CI.AtoD.allf20.lsdv.pearl <- pmax(0, pmin(1, CI.AtoD.allf20.lsdv.pearl))
CI.AtoD.new.lsdv.pearl <- pmax(0, pmin(1, CI.AtoD.new.lsdv.pearl))
CI.AtoD.newf10.lsdv.pearl <- pmax(0, pmin(1, CI.AtoD.newf10.lsdv.pearl))
CI.AtoD.newf20.lsdv.pearl <- pmax(0, pmin(1, CI.AtoD.newf20.lsdv.pearl))
CI.AtoD.old.lsdv.pearl <- pmax(0, pmin(1, CI.AtoD.old.lsdv.pearl))
CI.AtoD.oldf10.lsdv.pearl <- pmax(0, pmin(1, CI.AtoD.oldf10.lsdv.pearl))
CI.AtoD.oldf20.lsdv.pearl <- pmax(0, pmin(1, CI.AtoD.oldf20.lsdv.pearl))

CI.AtoD.allf10.lsdv.exog <- IMCI(l=resdem.AtoD.allf10.lsdv.exog[1], u=resdem.AtoD.allf10.lsdv.exog[2],
    var.l=var(resdem.AtoD.allf10.lsdv[,2], na.rm=T), var.u=var(resdem.AtoD.allf10.lsdv[,3], na.rm=T))
CI.AtoD.allf20.lsdv.exog <- IMCI(l=resdem.AtoD.allf20.lsdv.exog[1], u=resdem.AtoD.allf20.lsdv.exog[2],
    var.l=var(resdem.AtoD.allf20.lsdv[,2], na.rm=T), var.u=var(resdem.AtoD.allf20.lsdv[,3], na.rm=T))
CI.AtoD.new.lsdv.exog <- IMCI(l=resdem.AtoD.new.lsdv.exog[1], u=resdem.AtoD.new.lsdv.exog[2],
    var.l=var(resdem.AtoD.new.lsdv[,2], na.rm=T), var.u=var(resdem.AtoD.new.lsdv[,3], na.rm=T))
CI.AtoD.newf10.lsdv.exog <- IMCI(l=resdem.AtoD.newf10.lsdv.exog[1], u=resdem.AtoD.newf10.lsdv.exog[2],
    var.l=var(resdem.AtoD.newf10.lsdv[,2], na.rm=T), var.u=var(resdem.AtoD.newf10.lsdv[,3], na.rm=T))
CI.AtoD.newf20.lsdv.exog <- IMCI(l=resdem.AtoD.newf20.lsdv.exog[1], u=resdem.AtoD.newf20.lsdv.exog[2],
    var.l=var(resdem.AtoD.newf20.lsdv[,2], na.rm=T), var.u=var(resdem.AtoD.newf20.lsdv[,3], na.rm=T))
CI.AtoD.old.lsdv.exog <- IMCI(l=resdem.AtoD.old.lsdv.exog[1], u=resdem.AtoD.old.lsdv.exog[2],
    var.l=var(resdem.AtoD.old.lsdv[,2], na.rm=T), var.u=var(resdem.AtoD.old.lsdv[,3], na.rm=T))
CI.AtoD.oldf10.lsdv.exog <- IMCI(l=resdem.AtoD.oldf10.lsdv.exog[1], u=resdem.AtoD.oldf10.lsdv.exog[2],
    var.l=var(resdem.AtoD.oldf10.lsdv[,2], na.rm=T), var.u=var(resdem.AtoD.oldf10.lsdv[,3], na.rm=T))
CI.AtoD.oldf20.lsdv.exog <- IMCI(l=resdem.AtoD.oldf20.lsdv.exog[1], u=resdem.AtoD.oldf20.lsdv.exog[2],
    var.l=var(resdem.AtoD.oldf20.lsdv[,2], na.rm=T), var.u=var(resdem.AtoD.oldf20.lsdv[,3], na.rm=T))

CI.AtoD.allf10.lsdv.exog <- pmax(0, pmin(1, CI.AtoD.allf10.lsdv.exog$ci))
CI.AtoD.allf20.lsdv.exog <- pmax(0, pmin(1, CI.AtoD.allf20.lsdv.exog$ci))
CI.AtoD.new.lsdv.exog <- pmax(0, pmin(1, CI.AtoD.new.lsdv.exog$ci))
CI.AtoD.newf10.lsdv.exog <- pmax(0, pmin(1, CI.AtoD.newf10.lsdv.exog$ci))
CI.AtoD.newf20.lsdv.exog <- pmax(0, pmin(1, CI.AtoD.newf20.lsdv.exog$ci))
CI.AtoD.old.lsdv.exog <- pmax(0, pmin(1, CI.AtoD.old.lsdv.exog$ci))
CI.AtoD.oldf10.lsdv.exog <- pmax(0, pmin(1, CI.AtoD.oldf10.lsdv.exog$ci))
CI.AtoD.oldf20.lsdv.exog <- pmax(0, pmin(1, CI.AtoD.oldf20.lsdv.exog$ci))

################################################
# Figures
################################################

pdf("fig1-left.pdf", width=4, height=6)
plot(c(1,3), c(0,1), type="n", xlim=c(0.5,3.5), ylim=c(0,1), xaxt="n", xlab="",
    main = "Assassination and Regime Change", ylab="Probability of Causal Attribution")
abline(h=0, lty="dotted")
abline(h=0.2, lty="dotted")
abline(h=0.4, lty="dotted")
abline(h=0.6, lty="dotted")
abline(h=0.8, lty="dotted")
abline(h=1, lty="dotted")
mtext("All Changes", side=1, at=1, line=1, cex=.8)
mtext("Autocracy", side=1, at=2, line=1, cex=.8)
mtext("to", side=1, at=2, line=2, cex=.8)
mtext("Democracy", side=1, at=2, line=3, cex=.8)
mtext("Democracy", side=1, at=3, line=1, cex=.8)
mtext("to", side=1, at=3, line=2, cex=.8)
mtext("Autocracy", side=1, at=3, line=3, cex=.8)

# Point Estimates under Monotonicity
points(.9, resdem.all.lsdv.pearl, pch=16)
points(1.9, resdem.AtoD.lsdv.pearl, pch=16)
points(2.9, resdem.DtoA.lsdv.pearl, pch=16)

# Bounds under Exogeneity
lines(c(1.1,1.1), resdem.all.lsdv.exog, lwd=4)
lines(c(2.1,2.1), resdem.AtoD.lsdv.exog, lwd=4)
lines(c(3.1,3.1), resdem.DtoA.lsdv.exog, lwd=4)

# Confidence Intervals
lines(c(.9,.9), CI.all.lsdv.pearl)
lines(c(1.9,1.9), CI.AtoD.lsdv.pearl)
lines(c(2.9,2.9), CI.DtoA.lsdv.pearl)
points(c(.9,.9), CI.all.lsdv.pearl, pch="_", cex=2)
points(c(1.9,1.9), CI.AtoD.lsdv.pearl, pch="_", cex=2)
points(c(2.9,2.9), CI.DtoA.lsdv.pearl, pch="_", cex=2)

lines(c(1.1,1.1), CI.all.lsdv.exog)
lines(c(2.1,2.1), CI.AtoD.lsdv.exog)
lines(c(3.1,3.1), CI.DtoA.lsdv.exog)
points(c(1.1,1.1), CI.all.lsdv.exog, pch="_", cex=2)
points(c(2.1,2.1), CI.AtoD.lsdv.exog, pch="_", cex=2)
points(c(3.1,3.1), CI.DtoA.lsdv.exog, pch="_", cex=2)

dev.off()

pdf("fig1-right.pdf", width=6, height=6)
par(mfrow=c(3,3), oma=c(1,5,5,1), mar=c(.2,3,2.5,1))

plot(c(0.6,0.6), resdem.AtoD.lsdv.exog, type="l", lwd=4, xaxt="n",
    xlab = "", ylab = "", xlim = c(0,1), ylim = c(0,1))
points(0.4, resdem.AtoD.lsdv.pearl, pch=16, cex = 1.5)
lines(c(0.6,0.6), CI.AtoD.lsdv.exog)
points(c(0.6,0.6), CI.AtoD.lsdv.exog, pch="_", cex=2)
lines(c(0.4,0.4), CI.AtoD.lsdv.pearl)
points(c(0.4,0.4), CI.AtoD.lsdv.pearl, pch="_", cex=2)
abline(h=0, lty="dotted")
abline(h=0.2, lty="dotted")
abline(h=0.4, lty="dotted")
abline(h=0.6, lty="dotted")
abline(h=0.8, lty="dotted")
abline(h=1, lty="dotted")
mtext("All", side = 3, line = 1.3)
mtext("1 Year Out", side = 2, line = 3)

plot(c(0.6,0.6), resdem.AtoD.new.lsdv.exog, type="l", lwd=4, xaxt="n",
    xlab = "", ylab = "", xlim = c(0,1), ylim = c(0,1))
points(0.4, resdem.AtoD.new.lsdv.pearl, pch=16, cex = 1.5)
lines(c(0.6,0.6), CI.AtoD.new.lsdv.exog)
points(c(0.6,0.6), CI.AtoD.new.lsdv.exog, pch="_", cex=2)
lines(c(0.4,0.4), CI.AtoD.new.lsdv.pearl)
points(c(0.4,0.4), CI.AtoD.new.lsdv.pearl, pch="_", cex=2)
abline(h=0, lty="dotted")
abline(h=0.2, lty="dotted")
abline(h=0.4, lty="dotted")
abline(h=0.6, lty="dotted")
abline(h=0.8, lty="dotted")
abline(h=1, lty="dotted")
mtext(expression(paste("Tenure " <= 10)), side = 3, line = 1.3)


plot(c(0.6,0.6), resdem.AtoD.old.lsdv.exog, type="l", lwd=4, xaxt="n",
    xlab = "", ylab = "", xlim = c(0,1), ylim = c(0,1))
points(0.4, resdem.AtoD.old.lsdv.pearl, pch=16, cex = 1.5)
lines(c(0.6,0.6), CI.AtoD.old.lsdv.exog)
points(c(0.6,0.6), CI.AtoD.old.lsdv.exog, pch="_", cex=2)
lines(c(0.4,0.4), CI.AtoD.old.lsdv.pearl)
points(c(0.4,0.4), CI.AtoD.old.lsdv.pearl, pch="_", cex=2)
abline(h=0, lty="dotted")
abline(h=0.2, lty="dotted")
abline(h=0.4, lty="dotted")
abline(h=0.6, lty="dotted")
abline(h=0.8, lty="dotted")
abline(h=1, lty="dotted")
mtext(expression(paste("Tenure " > 10)), side = 3, line = 1.3)

plot(c(0.6,0.6), resdem.AtoD.allf10.lsdv.exog, type="l", lwd=4, xaxt="n",
    xlab = "", ylab = "", xlim = c(0,1), ylim = c(0,1))
points(0.4, resdem.AtoD.allf10.lsdv.pearl, pch=16, cex = 1.5)
lines(c(0.6,0.6), CI.AtoD.allf10.lsdv.exog)
points(c(0.6,0.6), CI.AtoD.allf10.lsdv.exog, pch="_", cex=2)
lines(c(0.4,0.4), CI.AtoD.allf10.lsdv.pearl)
points(c(0.4,0.4), CI.AtoD.allf10.lsdv.pearl, pch="_", cex=2)
abline(h=0, lty="dotted")
abline(h=0.2, lty="dotted")
abline(h=0.4, lty="dotted")
abline(h=0.6, lty="dotted")
abline(h=0.8, lty="dotted")
abline(h=1, lty="dotted")
mtext("10 Years Out", side = 2, line = 3)

plot(c(0.6,0.6), resdem.AtoD.newf10.lsdv.exog, type="l", lwd=4, xaxt="n",
    xlab = "", ylab = "", xlim = c(0,1), ylim = c(0,1))
points(0.4, resdem.AtoD.newf10.lsdv.pearl, pch=16, cex = 1.5)
lines(c(0.6,0.6), CI.AtoD.newf10.lsdv.exog)
points(c(0.6,0.6), CI.AtoD.newf10.lsdv.exog, pch="_", cex=2)
lines(c(0.4,0.4), CI.AtoD.newf10.lsdv.pearl)
points(c(0.4,0.4), CI.AtoD.newf10.lsdv.pearl, pch="_", cex=2)
abline(h=0, lty="dotted")
abline(h=0.2, lty="dotted")
abline(h=0.4, lty="dotted")
abline(h=0.6, lty="dotted")
abline(h=0.8, lty="dotted")
abline(h=1, lty="dotted")

plot(c(0.6,0.6), resdem.AtoD.oldf10.lsdv.exog, type="l", lwd=4, xaxt="n",
    xlab = "", ylab = "", xlim = c(0,1), ylim = c(0,1))
points(0.4, resdem.AtoD.oldf10.lsdv.pearl, pch=16, cex = 1.5)
lines(c(0.6,0.6), CI.AtoD.oldf10.lsdv.exog)
points(c(0.6,0.6), CI.AtoD.oldf10.lsdv.exog, pch="_", cex=2)
lines(c(0.4,0.4), CI.AtoD.oldf10.lsdv.pearl)
points(c(0.4,0.4), CI.AtoD.oldf10.lsdv.pearl, pch="_", cex=2)
abline(h=0, lty="dotted")
abline(h=0.2, lty="dotted")
abline(h=0.4, lty="dotted")
abline(h=0.6, lty="dotted")
abline(h=0.8, lty="dotted")
abline(h=1, lty="dotted")

plot(c(0.6,0.6), resdem.AtoD.allf20.lsdv.exog, type="l", lwd=4, xaxt="n",
    xlab = "", ylab = "", xlim = c(0,1), ylim = c(0,1))
points(0.4, resdem.AtoD.allf20.lsdv.pearl, pch=16, cex = 1.5)
lines(c(0.6,0.6), CI.AtoD.allf20.lsdv.exog)
points(c(0.6,0.6), CI.AtoD.allf20.lsdv.exog, pch="_", cex=2)
lines(c(0.4,0.4), CI.AtoD.allf20.lsdv.pearl)
points(c(0.4,0.4), CI.AtoD.allf20.lsdv.pearl, pch="_", cex=2)
abline(h=0, lty="dotted")
abline(h=0.2, lty="dotted")
abline(h=0.4, lty="dotted")
abline(h=0.6, lty="dotted")
abline(h=0.8, lty="dotted")
abline(h=1, lty="dotted")
mtext("20 Years Out", side = 2, line = 3)

plot(c(0.6,0.6), resdem.AtoD.newf20.lsdv.exog, type="l", lwd=4, xaxt="n",
    xlab = "", ylab = "", xlim = c(0,1), ylim = c(0,1))
points(0.4, resdem.AtoD.newf20.lsdv.pearl, pch=16, cex = 1.5)
lines(c(0.6,0.6), CI.AtoD.newf20.lsdv.exog)
points(c(0.6,0.6), CI.AtoD.newf20.lsdv.exog, pch="_", cex=2)
lines(c(0.4,0.4), CI.AtoD.newf20.lsdv.pearl)
points(c(0.4,0.4), CI.AtoD.newf20.lsdv.pearl, pch="_", cex=2)
abline(h=0, lty="dotted")
abline(h=0.2, lty="dotted")
abline(h=0.4, lty="dotted")
abline(h=0.6, lty="dotted")
abline(h=0.8, lty="dotted")
abline(h=1, lty="dotted")

plot(c(0.6,0.6), resdem.AtoD.oldf20.lsdv.exog, type="l", lwd=4, xaxt="n",
    xlab = "", ylab = "", xlim = c(0,1), ylim = c(0,1))
points(0.4, resdem.AtoD.oldf20.lsdv.pearl, pch=16, cex = 1.5)
lines(c(0.6,0.6), CI.AtoD.oldf20.lsdv.exog)
points(c(0.6,0.6), CI.AtoD.oldf20.lsdv.exog, pch="_", cex=2)
lines(c(0.4,0.4), CI.AtoD.oldf20.lsdv.pearl)
points(c(0.4,0.4), CI.AtoD.oldf20.lsdv.pearl, pch="_", cex=2)
abline(h=0, lty="dotted")
abline(h=0.2, lty="dotted")
abline(h=0.4, lty="dotted")
abline(h=0.6, lty="dotted")
abline(h=0.8, lty="dotted")
abline(h=1, lty="dotted")

title(main = "Assassination and Democratic Transition",
    outer = T, cex.main = 2, line = 3)
mtext("Tenure of Leader", side = 3, outer = T, cex = 1.25, line = .8)
mtext("Time after Transition", side = 2, outer = T, line = 3, cex = 1.25)

dev.off()

